/****************************************************************************
 * Copyright (C) 2009-2010 SciTouch LLC
 * 
 * This file is part of Indigo toolkit.
 * 
 * This file may be distributed and/or modified under the terms of the
 * GNU General Public License version 3 as published by the Free Software
 * Foundation and appearing in the file LICENSE.GPL included in the
 * packaging of this file.
 * 
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 ***************************************************************************/

#include "molecule/molecule.h"
#include "molecule/molecule_arom.h"
#include "base_cpp/gray_codes.h"
#include "base_c/bitarray.h"

//
// AromatizerBase
// 
AromatizerBase::AromatizerBase (const Molecule &molecule) : _molecule(molecule),
   TL_CP_GET(_bonds_arom),
   TL_CP_GET(_bonds_arom_count),
   TL_CP_GET(_unsure_cycles)
{
   _bonds_arom.resize(bitGetSize(molecule.edgeEnd()));
   _bonds_arom_count.resize(molecule.edgeEnd());
   reset();
}

AromatizerBase::~AromatizerBase ()
{
}

bool AromatizerBase::_checkVertex (int v_idx)
{
   return true;
}

bool AromatizerBase::_isCycleAromatic (const int *cycle, int cycle_len)
{
   return true;
}

bool AromatizerBase::_checkDoubleBonds (const int *cycle, int cycle_len)
{
   int j;

   for (j = 0; j < cycle_len; j++)
   {
      int v_left_idx = cycle[j];
      int v_center_idx = cycle[(j + 1) % cycle_len];
      int v_right_idx = cycle[(j + 2) % cycle_len];

      const Vertex &vertex = _molecule.getVertex(v_center_idx);
      int i;

      int internalDoubleBondCount = 0;
      for (i = vertex.neiBegin(); i != vertex.neiEnd(); i = vertex.neiNext(i))
      {
         int nei_idx = vertex.neiVertex(i);

         int e_idx = vertex.neiEdge(i);
         int type = _molecule.getBond(e_idx).type;
         if (type == BOND_DOUBLE && !isBondAromatic(e_idx)) {
            if (nei_idx != v_left_idx && nei_idx != v_right_idx)
               return false;
            else if (nei_idx == v_left_idx)
               internalDoubleBondCount++;
            else if (nei_idx == v_left_idx)
               internalDoubleBondCount++;
         }
      }
      if (internalDoubleBondCount >= 2)
         return false;
   }
   return true;
}

void AromatizerBase::_aromatizeCycle (const int *cycle, int cycle_len)
{
   for (int i = 0; i < cycle_len; i++)
   {
      int a = cycle[i], b = cycle[(i + 1) % cycle_len];
      int e_idx = _molecule.findEdgeIndex(a, b);
      _bonds_arom_count[e_idx]++;
      bitSetBit(_bonds_arom.ptr(), e_idx, 1);
   }
   _handleAromaticCycle(cycle, cycle_len);
}

void AromatizerBase::_handleCycle (const Array<int> &path)
{
   // Check Huckel's rule
   if (!_isCycleAromatic(path.ptr(), path.size()))
      return;
   addAromaticCycle(-1, path.ptr(), path.size());
}

bool AromatizerBase::handleUnsureCycles ()
{
   bool changed = true;
   bool is_all_aromatic = true;

   while (changed)
   {
      changed = false;

      for (int i = 0; i < _unsure_cycles.size(); i++)
      {
         if (_unsure_cycles[i].isEmpty)
            continue;

         if (_checkDoubleBonds(_unsure_cycles[i].cycle, _unsure_cycles[i].length))
         {
            _aromatizeCycle(_unsure_cycles[i].cycle, _unsure_cycles[i].length);
            _unsure_cycles[i].isEmpty = true;
            changed = true;
         }
         else 
            is_all_aromatic = false;
      }
   }
   return is_all_aromatic;
}

bool AromatizerBase::_cb_check_vertex (const Graph &graph, int v_idx, void *context)
{
   AromatizerBase *arom = (AromatizerBase *)context;
   return arom->_checkVertex(v_idx);
}

bool AromatizerBase::_cb_handle_cycle (const Graph &graph, const Array<int> &vertices, const Array<int> &edges, void *context)
{
   AromatizerBase *arom = (AromatizerBase *)context;
   arom->_handleCycle(vertices);
   return true;
}

void AromatizerBase::aromatize ()
{
   CycleEnumerator cycle_enumerator(_molecule);

   cycle_enumerator.cb_check_vertex = _cb_check_vertex;
   cycle_enumerator.cb_handle_cycle = _cb_handle_cycle;
   cycle_enumerator.max_length = MAX_CYCLE_LEN;
   cycle_enumerator.context = this;
   cycle_enumerator.process();

   handleUnsureCycles();
}

bool AromatizerBase::isBondAromatic (int e_idx)
{
   return _bonds_arom_count[e_idx] != 0;
}

const byte* AromatizerBase::isBondAromaticArray (void)
{
   return _bonds_arom.ptr();
}

void AromatizerBase::addAromaticCycle (int id, const int *cycle, int cycle_len)
{
   if (!_checkDoubleBonds(cycle, cycle_len))
   {
      int empty_idx;
      if (_unsureCyclesCount == _unsure_cycles.size())
      {
         empty_idx = _unsure_cycles.size();
         _unsure_cycles.push();
      }
      else
      {
         // Find first empty space
         empty_idx = -1;
         for (int i = 0; i < _unsure_cycles.size(); i++)
            if (_unsure_cycles[i].isEmpty)
            {
               empty_idx = i;
               break;
            }
         if (empty_idx == -1)
            throw Exception("AromatizerBase::addAromaticCycle: internal logic error");
      }

      CycleDef &cycleDef = _unsure_cycles[empty_idx];
      cycleDef.id = id;
      cycleDef.isEmpty = false;
      cycleDef.length = cycle_len;

      memcpy(cycleDef.cycle, cycle, cycle_len * sizeof(int));

      _unsureCyclesCount++;
   }
   else
      _aromatizeCycle(cycle, cycle_len);
}

void AromatizerBase::removeAromaticCycle (int id, const int *cycle, int cycle_len)
{
   // Find id in unsure cycles
   for (int i = 0; i < _unsure_cycles.size(); i++)
      if (!_unsure_cycles[i].isEmpty && _unsure_cycles[i].id == id)
      {
         _unsure_cycles[i].isEmpty = true;
         return;
      }

   // Just decrease bond marks
   for (int i = 0; i < cycle_len; i++)
   {
      int a = cycle[i], b = cycle[(i + 1) % cycle_len];
      int e_idx = _molecule.findEdgeIndex(a, b);
      _bonds_arom_count[e_idx]--;
      if (_bonds_arom_count[e_idx] == 0)
         bitSetBit(_bonds_arom.ptr(), e_idx, 1);
   }
}

void AromatizerBase::_handleAromaticCycle (const int *cycle, int cycle_len)
{
}

void AromatizerBase::reset (void)
{
   _unsure_cycles.clear();
   _bonds_arom.zerofill();
   _bonds_arom_count.zerofill();
   _cyclesHandled = 0;
   _unsureCyclesCount = 0;
}

void AromatizerBase::setBondAromaticCount (int e_idx, int count)
{
   _bonds_arom_count[e_idx] = count;
   bitSetBit(_bonds_arom.ptr(), e_idx, count != 0);
}

//
// MoleculeAromatizer
// 

MoleculeAromatizer::MoleculeAromatizer (const Molecule &molecule) : AromatizerBase(molecule),
   TL_CP_GET(_pi_labels)
{
   _pi_labels.clear_resize(molecule.vertexEnd());
}

int MoleculeAromatizer::_getPiLabel (int v_idx)
{
   const Vertex &vertex = _molecule.getVertex(v_idx);
   int lonepairs = 0;
   const Atom &atom = _molecule.getAtom(v_idx);
   int radical = atom.radical;
   int i;

   if (_molecule.isPseudoAtom(v_idx))
      return -1;

   if (!Element::canBeAromatic(atom.label))
      return -1;
   
   if (radical > 0)
      return 1;

   for (i = vertex.neiBegin(); i != vertex.neiEnd(); i = vertex.neiNext(i))
   {
      int type = _molecule.getBond(vertex.neiEdge(i)).type;

      if (type == BOND_DOUBLE/* || type == BOND_TRIPLE*/)
         return 1;
   }

   if (_molecule.getVacantPiOrbitals(v_idx, &lonepairs) > 0)
      return 0;

   if (lonepairs > 0)
      return 2;

   return -1;
}

void MoleculeAromatizer::precalculatePiLabels ()
{
   for (int v_idx = _molecule.vertexBegin(); 
         v_idx < _molecule.vertexEnd(); 
         v_idx = _molecule.vertexNext(v_idx))
      _pi_labels[v_idx] = _getPiLabel(v_idx);
}

bool MoleculeAromatizer::_checkVertex (int v_idx)
{
   return _pi_labels[v_idx] != -1;
}

bool MoleculeAromatizer::_isCycleAromatic (const int *cycle, int cycle_len)
{
   int count = 0;
   // Check Huckel's rule
   for (int i = 0; i < cycle_len; i++)
      count += _pi_labels[cycle[i]];

   if (((count - 2) % 4) != 0) 
      return false;
   return true;
}

void MoleculeAromatizer::aromatizeBonds (Molecule &mol)
{
   MoleculeAromatizer aromatizer(mol);

   aromatizer.precalculatePiLabels();
   aromatizer.aromatize();

   for (int e_idx = mol.edgeBegin(); e_idx < mol.edgeEnd(); e_idx = mol.edgeNext(e_idx))
      if (aromatizer.isBondAromatic(e_idx))
         mol.setBondAromatic(e_idx);
}

//
// QueryMoleculeAromatizer
//

QueryMoleculeAromatizer::QueryMoleculeAromatizer (const Molecule &molecule) : AromatizerBase(molecule), 
   TL_CP_GET(_pi_labels),
   TL_CP_GET(_aromatic_cycles),
   TL_CP_GET(_result_molecule)
{
   _pi_labels.clear_resize(molecule.vertexEnd());
   _aromatic_cycles.clear();
   _aromatic_cycles.reserve(100);
   _mode = FUZZY;
   _collecting = false;
}

void QueryMoleculeAromatizer::precalculatePiLabels ()
{
   for (int v_idx = _molecule.vertexBegin(); 
            v_idx < _molecule.vertexEnd(); 
            v_idx = _molecule.vertexNext(v_idx))
      _pi_labels[v_idx] = _getPiLabel(v_idx);
}

bool QueryMoleculeAromatizer::_checkVertex (int v_idx)
{
   return _pi_labels[v_idx].min != -1;
}

bool QueryMoleculeAromatizer::_isCycleAromatic (const int *cycle, int cycle_len)
{
   // Single/double bond can't be aromatic and Check if cycle wasn't aromatic
   bool all_aromatic = true;
   for (int i = 0; i < cycle_len; i++)
   {
      int a = cycle[i], b = cycle[(i + 1) % cycle_len];
      int e_idx = _molecule.findEdgeIndex(a, b);
      const QueryBond &query_bond = _molecule.getQueryBond(e_idx);
      if (query_bond.type == QUERY_BOND_SINGLE_OR_DOUBLE)
         return false;
      const Bond &bond = _molecule.getBond(e_idx);
      if (query_bond.type != 0 || bond.type != BOND_AROMATIC)
         all_aromatic = false;
   }
   if (all_aromatic)
      return false;

   PiValue cycleSum(0, 0, true);
   // Check Huckel's rule
   for (int i = 0; i < cycle_len; i++)
   {
      PiValue &cur = _pi_labels[cycle[i]];
      if (cur.max != -1)
      {
         cycleSum.max += cur.max;
         cycleSum.min += cur.min;
      }
      else
      {
         cycleSum.max += cur.max;
         cycleSum.min += cur.min;
      }
      if (cur.exact != -1 && cycleSum.exact != -1)
         cycleSum.exact += cur.exact;
      else
         cycleSum.exact = -1;
   }

   // Check Huckel's rule
   if (_mode == EXACT)
   {
      int sum = 0;
      for (int i = 0; i < cycle_len; i++)
      {
         int cur = _pi_labels[cycle[i]].exact;
         if (cur != -1)
            sum += cur;
         else
         {
            sum = -1;
            break;
         }
      }
      // Check if cycle have pi-lables sum 4n+2 for drawn query
      if (sum == -1 || sum % 4 != 2)
         return false;
      return true;
   }

   //
   // Fuzzy mode: check if circle can have 4n-2 value
   //

   if (cycleSum.max - cycleSum.min > 3)
      return true;

   int cycle_min = 0;
   int cycle_max = 0;
   for (int i = 0; i < cycle_len; i++)
   {
      PiValue &cur = _pi_labels[cycle[i]];
      if (cur.min != -1)
      {
         cycle_min += cur.min;
         cycle_max += cur.max;
      }
      else
      {
         cycle_min = -1;
         break;
      }
   }

   if (cycle_min == -1)
      return false;

   int residue_min = (cycle_min + 2) % 4;
   int residue_max = (cycle_max + 2) % 4;

   if (residue_min == 0 || residue_min > residue_max) 
      return true;
   return false;
}

QueryMoleculeAromatizer::PiValue QueryMoleculeAromatizer::_getPiLabel (int v_idx)
{
   const Vertex &vertex = _molecule.getVertex(v_idx);
   const Atom &atom = _molecule.getAtom(v_idx);
   const QueryAtom &query_atom = _molecule.getQueryAtom(v_idx);
   bool 
      have_double_bond = false, 
      have_query_bond = false;

   if (atom.radical > 0)
      return PiValue(1, 1, true);

   // Check double bonds
   int min_conn = vertex.degree();
   for (int i = vertex.neiBegin(); i != vertex.neiEnd(); i = vertex.neiNext(i))
   {
      const Bond &bond = _molecule.getBond(vertex.neiEdge(i));
      const QueryBond &queryBond = _molecule.getQueryBond(vertex.neiEdge(i));

      if (queryBond.type != 0)
         have_query_bond = true;
      else
         switch (bond.type)
         {
         case BOND_DOUBLE:
            have_double_bond = true;
            min_conn++;
            break;
         case BOND_TRIPLE:
            // Triple bonds can't be attached to atoms in aromatic rings
            return PiValue(-1, -1, true);
         case BOND_AROMATIC:
            have_query_bond = true;
            break;
         }
   }

   if (query_atom.type != 0)
   {
      if (query_atom.type == QUERY_ATOM_RGROUP)
      {
         if (vertex.degree() == 1)
         {
            // R-Group with single attachment point
            if (have_double_bond)
               return PiValue(1, 1);
            return PiValue(0, 2);
         }
         else
         {
            // R-Group with two attachment points
            // Here can be chain of atoms with different summary pi-labels
            return PiValue(0, 4);
         }
      }
      else
      {
         if (have_double_bond)
         {
            int exact_aromatic = -1;
            if (query_atom.type == QUERY_ATOM_LIST)
            {
               exact_aromatic = 1;
               for (int i = 0; i < query_atom.list_size; i++)
                  if (!Element::canBeAromatic(query_atom.list[i]))
                  {
                     exact_aromatic = -1;
                     break;
                  }
            }
            return PiValue(1, 1, exact_aromatic);
         }
         return PiValue(0, 2);
      }
   }
   int max_conn = Element::getMaximumConnectivity(atom.label, 
      atom.charge, atom.radical, false);

   if (max_conn < min_conn && !have_double_bond && !have_query_bond)
   {
      // Such atom can't in aromatic ring (P with 5 valence for example)
      return PiValue(-1, -1, -1);
   }

   if (have_double_bond)
   {
      // Atom can be in aromatic ring with pi label = 1
      if (!have_query_bond)
      {
         // If drawn query molecule is aromatic then atom pi label = 1
         return PiValue(1, 1, 1);
      }
      else
      {
         // Atom can be or not in aromatic ring but pi label = 1
         return PiValue(1, 1, -1);
      }
   }

   int val, hyd;
   Element::calcValence(atom.label, 
      atom.charge, atom.radical, min_conn, val, hyd, true);

   int pi_label_without_double_bound = -1;
   int lonepairs;
   if (_molecule.getVacantPiOrbitals(v_idx, min_conn + hyd, &lonepairs) > 0)
      pi_label_without_double_bound = 0;
   else if (lonepairs > 0)
      pi_label_without_double_bound = 2;

   bool can_have_double_bond = have_query_bond;
   if (!can_have_double_bond)
   {
      int pi = pi_label_without_double_bound;
      return PiValue(pi, pi, pi);
   }

   // Atom have query bond attached
   if (pi_label_without_double_bound == -1)
      return PiValue(1, 1, -1);
   if (pi_label_without_double_bound == 0)
      return PiValue(0, 1, -1);
   else /* (pi_label_without_double_bound == 2) */
      return PiValue(1, 2, -1);
}

bool QueryMoleculeAromatizer::enumerateMolecules ()
{
   setMode(QueryMoleculeAromatizer::FUZZY);

   precalculatePiLabels();
   // Collect all aromatic cycles
   _collecting = true;
   aromatize();
   _collecting = false;
   reset();

   // Check estimation number of results molecules
   if (cb_check_cycles_count && !cb_check_cycles_count(_aromatic_cycles.size(), context))
      return false;

   Molecule &qmol = (Molecule &)_molecule;

   // Do nothing
   if (_aromatic_cycles.size() == 0)
   {
      if (cb_handle_molecule_aromatization)
         cb_handle_molecule_aromatization(qmol, context);
      return true;
   }

   // Set aromatic bonds aromatic count to infinity
   for (int e_idx = qmol.edgeBegin(); 
      e_idx < qmol.edgeEnd(); 
      e_idx = qmol.edgeNext(e_idx))
   {
      const QueryBond &query_bond = qmol.getQueryBond(e_idx);
      const Bond &bond = qmol.getBond(e_idx);
      if (query_bond.type == 0 && bond.type == BOND_AROMATIC)
         setBondAromaticCount(e_idx, 10000);
   }

   // LIFO
   enum
   {
      RECORDS_COUNT = 50
   };
   
   int proceededCount = 0;
   int currentIndex = 0;
   int recordSize = bitGetSize(qmol.edgeEnd());
   QS_DEF(Array<byte>, records);
   records.resize(RECORDS_COUNT * recordSize);

   // Enumerate all configurations
   GrayCodesEnumerator grayCodes(_aromatic_cycles.size(), true);
   do
   {
      if (handleUnsureCycles())
      {
         // Check if this state is the same as previous
         int end = __min(proceededCount, RECORDS_COUNT);
         const byte *aromState = isBondAromaticArray();
         bool same = false;
         for (int i = 0; i < end; i++)
            if (bitTestEquality(aromState, &records[i * recordSize], qmol.edgeEnd()))
            {
               same = true;
               break;
            }
         if (!same)
         {
            // Add record
            proceededCount++;
            memcpy(&records[currentIndex * recordSize], aromState, recordSize);
            currentIndex = (currentIndex + 1) % RECORDS_COUNT;

            // handle this aromatic configuration
            QS_DEF(Array<int>, inv_mapping);
            _result_molecule.cloneMolecule(qmol, NULL, &inv_mapping);
            for (int e_idx = qmol.edgeBegin(); 
               e_idx < qmol.edgeEnd(); 
               e_idx = qmol.edgeNext(e_idx))
            {
               int e_idx_mapped = Graph::findMappedEdge(qmol, 
                  _result_molecule, e_idx, inv_mapping.ptr());

               if (isBondAromatic(e_idx))
                  _result_molecule.setQueryBondAromatic(e_idx_mapped);
            }
            if (cb_handle_molecule_aromatization)
               if (!cb_handle_molecule_aromatization(_result_molecule, context))
                  return true;
         }
     }

      grayCodes.next();

      if (!grayCodes.isDone())
      {
         int c_idx = grayCodes.getBitChangeIndex();
         const byte *code = grayCodes.getCode();
         int need_add = bitGetBit(code, c_idx);
         if (need_add)
            addAromaticCycle(c_idx, _aromatic_cycles[c_idx].cycle, _aromatic_cycles[c_idx].length);
         else
            removeAromaticCycle(c_idx, _aromatic_cycles[c_idx].cycle, _aromatic_cycles[c_idx].length);
      }
   } while (!grayCodes.isDone());

   return true;
}

void QueryMoleculeAromatizer::_handleAromaticCycle (const int *cycle, int cycle_len)
{
   if (!_collecting)
      return;
   // Add cycle to storage      
   _aromatic_cycles.push();
   CycleDef &def = _aromatic_cycles[_aromatic_cycles.size() - 1];
   def.id = _aromatic_cycles.size() - 1;
   def.isEmpty = false;
   def.length = cycle_len;
   memcpy(def.cycle, cycle, cycle_len * sizeof(int));

   AromatizerBase::_handleAromaticCycle(cycle, cycle_len);
}

void QueryMoleculeAromatizer::setMode (int mode)
{
   _mode = mode;
}

void QueryMoleculeAromatizer::aromatizeBonds (Molecule &mol)
{
   // Mark edges that can be aromatic in some matching
   _aromatizeBondsFuzzy(mol);
   // Aromatize all aromatic cycles
   _aromatizeBondsExact(mol);

   if (mol.getRGroups() != 0)
   {
      MoleculeRGroups &rgroups = *mol.getRGroups();
      int n_rgroups = rgroups.getRGroupCount();

      for (int i = 0; i < n_rgroups; i++)
      {
         RGroup &rgroup = mol.getRGroups()->getRGroup(i);
         if (rgroup.fragments.size() > 0)
         {
            for (int j = 0; j < rgroup.fragments.size(); j++)
            {
               Molecule &fragment = *rgroup.fragments[j];
               _aromatizeRGroupFragment(fragment);
            }
         }
      }
   }
}

void QueryMoleculeAromatizer::_aromatizeRGroupFragment (Molecule &fragment)
{
   // Add additional atom to attachment points
   int additional_atom = fragment.addAtom(0);
   fragment.getQueryAtom2(additional_atom).type = QUERY_ATOM_RGROUP;

   MoleculeRGroupFragment *rgroupFragment = 
      fragment.getRGroupFragment();

   // Connect it with attachment points
   int maxOrder = rgroupFragment->attachmentPointCount();
   for (int i = 0; i < maxOrder; i++) 
   {
      int pointIndex = 0;
      int point;
      while (true)
      {
         point = rgroupFragment->getAttachmentPoint(i, pointIndex);
         if (point == -1)
            break;

         if (fragment.findEdgeIndex(point, additional_atom) == -1)
         {
            int qb = fragment.addBond(point, additional_atom, 0);
            fragment.getQueryBond2(qb).type = QUERY_BOND_ANY;
         }

         pointIndex++;
      }
   }

   aromatizeBonds(fragment);

   QS_DEF(Array<int>, indices);
   indices.clear();
   indices.push(additional_atom);
   fragment.removeAtoms(indices);
}

// Some cycles with query features can be aromatized
void QueryMoleculeAromatizer::_aromatizeBondsExact (Molecule &qmol)
{
   QueryMoleculeAromatizer aromatizer(qmol);

   aromatizer.setMode(QueryMoleculeAromatizer::EXACT);
   aromatizer.precalculatePiLabels();
   aromatizer.aromatize();

   for (int e_idx = qmol.edgeBegin(); e_idx < qmol.edgeEnd(); e_idx = qmol.edgeNext(e_idx))
      if (aromatizer.isBondAromatic(e_idx))
         qmol.setQueryBondAromatic(e_idx);
}

void QueryMoleculeAromatizer::_aromatizeBondsFuzzy (Molecule &mol)
{
   QueryMoleculeAromatizer aromatizer(mol);

   aromatizer.setMode(QueryMoleculeAromatizer::FUZZY);
   aromatizer.precalculatePiLabels();
   aromatizer.aromatize();

   for (int e_idx = mol.edgeBegin(); e_idx < mol.edgeEnd(); e_idx = mol.edgeNext(e_idx))
      if (aromatizer.isBondAromatic(e_idx))
         mol.setQueryBondFuzzyAromatic(e_idx);
}

void MoleculeAromatizer::findAromaticAtoms (const Molecule &mol, Array<int> *atoms, Array<int> *bonds, bool allow_arom)
{
   Molecule clone;
   Array<int> mapping;

   if (atoms != 0)
   {
      atoms->clear_resize(mol.vertexEnd());
      atoms->zerofill();
   }

   if (bonds != 0)
   {
      bonds->clear_resize(mol.edgeEnd());
      bonds->zerofill();
   }

   clone.cloneMolecule(mol, &mapping, 0);

   clone.calcImplicitHydrogens(allow_arom);
   clone.findBondsInRings();

   MoleculeAromatizer::aromatizeBonds(clone);

   for (int i = clone.edgeBegin(); i != clone.edgeEnd(); i = clone.edgeNext(i))
   {
      if (clone.getBond(i).type == BOND_AROMATIC)
      {
         const Edge &edge = clone.getEdge(i);

         if (atoms != 0)
         {
            atoms->at(mapping[edge.beg]) = 1;
            atoms->at(mapping[edge.end]) = 1;
         }

         if (bonds != 0)
            bonds->at(mol.findEdgeIndex(mapping[edge.beg], mapping[edge.end])) = 1;

      }
   }
}
